{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       "@import url('http://fonts.googleapis.com/css?family=Source+Code+Pro');\n",
       "@import url('http://fonts.googleapis.com/css?family=Vollkorn');\n",
       "@import url('http://fonts.googleapis.com/css?family=Arimo');\n",
       "@import url('http://fonts.googleapis.com/css?family=Fira_sans');\n",
       "\n",
       "    div.cell{\n",
       "        width: 900px;\n",
       "        margin-left: 0% !important;\n",
       "        margin-right: auto;\n",
       "    }\n",
       "    div.text_cell code {\n",
       "        background: transparent;\n",
       "        color: #000000;\n",
       "        font-weight: 600;\n",
       "        font-size: 11pt;\n",
       "        font-style: bold;\n",
       "        font-family:  'Source Code Pro', Consolas, monocco, monospace;\n",
       "   }\n",
       "    h1 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "\t}\n",
       "\t\n",
       "    div.input_area {\n",
       "        background: #F6F6F9;\n",
       "        border: 1px solid #586e75;\n",
       "    }\n",
       "\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 30pt;\n",
       "        line-height: 100%;\n",
       "        color:#c76c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 1em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    } \n",
       "    h2 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "        text-align: left;\n",
       "    }\n",
       "    .text_cell_render h2 {\n",
       "        font-weight: 200;\n",
       "        font-size: 16pt;\n",
       "        font-style: italic;\n",
       "        line-height: 100%;\n",
       "        color:#c76c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 1.5em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    } \n",
       "    h3 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h3 {\n",
       "        font-weight: 200;\n",
       "        font-size: 14pt;\n",
       "        line-height: 100%;\n",
       "        color:#d77c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 2em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    }\n",
       "    h4 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h4 {\n",
       "        font-weight: 100;\n",
       "        font-size: 14pt;\n",
       "        color:#d77c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "        white-space: nowrap;\n",
       "    }\n",
       "    h5 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 200;\n",
       "        font-style: normal;\n",
       "        color: #1d3b84;\n",
       "        font-size: 16pt;\n",
       "        margin-bottom: 0em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "        white-space: nowrap;\n",
       "    }\n",
       "    div.text_cell_render{\n",
       "        font-family: 'Fira sans', verdana,arial,sans-serif;\n",
       "        line-height: 125%;\n",
       "        font-size: 115%;\n",
       "        text-align:justify;\n",
       "        text-justify:inter-word;\n",
       "    }\n",
       "    div.output_subarea.output_text.output_pyout {\n",
       "        overflow-x: auto;\n",
       "        overflow-y: scroll;\n",
       "        max-height: 50000px;\n",
       "    }\n",
       "    div.output_subarea.output_stream.output_stdout.output_text {\n",
       "        overflow-x: auto;\n",
       "        overflow-y: scroll;\n",
       "        max-height: 50000px;\n",
       "    }\n",
       "    div.output_wrapper{\n",
       "        margin-top:0.2em;\n",
       "        margin-bottom:0.2em;\n",
       "}\n",
       "\n",
       "    code{\n",
       "      font-size: 70%;\n",
       "    }\n",
       "    .rendered_html code{\n",
       "    background-color: transparent;\n",
       "    }\n",
       "    ul{\n",
       "        margin: 2em;\n",
       "    }\n",
       "    ul li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.5em; \n",
       "    }\n",
       "    ul li li{\n",
       "        padding-left: 0.2em; \n",
       "        margin-bottom: 0.2em; \n",
       "        margin-top: 0.2em; \n",
       "    }\n",
       "    ol{\n",
       "        margin: 2em;\n",
       "    }\n",
       "    ol li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.5em; \n",
       "    }\n",
       "    ul li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.2em; \n",
       "    }\n",
       "    a:link{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    a:visited{\n",
       "       font-weight: bold;\n",
       "       color: #1d3b84;\n",
       "    }\n",
       "    a:hover{\n",
       "       font-weight: bold;\n",
       "       color: #1d3b84;\n",
       "    }\n",
       "    a:focus{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    a:active{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    .rendered_html :link {\n",
       "       text-decoration: underline; \n",
       "    }\n",
       "    .rendered_html :hover {\n",
       "       text-decoration: none; \n",
       "    }\n",
       "    .rendered_html :visited {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .rendered_html :focus {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .rendered_html :active {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "    } \n",
       "    hr {\n",
       "      color: #f3f3f3;\n",
       "      background-color: #f3f3f3;\n",
       "      height: 1px;\n",
       "    }\n",
       "    blockquote{\n",
       "      display:block;\n",
       "      background: #fcfcfc;\n",
       "      border-left: 5px solid #c76c0c;\n",
       "      font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "      width:680px;\n",
       "      padding: 10px 10px 10px 10px;\n",
       "      text-align:justify;\n",
       "      text-justify:inter-word;\n",
       "      }\n",
       "      blockquote p {\n",
       "        margin-bottom: 0;\n",
       "        line-height: 125%;\n",
       "        font-size: 100%;\n",
       "      }\n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    scale:100,\n",
       "                        availableFonts: [],\n",
       "                        preferredFont:null,\n",
       "                        webFont: \"TeX\",\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# format the book\n",
    "%matplotlib inline\n",
    "import sys\n",
    "from __future__ import division, print_function\n",
    "import sys\n",
    "sys.path.insert(0,'../code')\n",
    "import book_format\n",
    "book_format.load_style('../code')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation\n",
    "\n",
    "In this chapter I describe my solution to a problem posed by a patient\n",
    "with a kidney tumor. I think the problem is important and relevant to\n",
    "patients with these tumors and doctors treating them.\n",
    "\n",
    "And I think the solution is interesting because, although it is a\n",
    "Bayesian approach to the problem, the use of Bayes’s theorem is\n",
    "implicit. I present the solution and my code; at the end of the chapter\n",
    "I will explain the Bayesian part.\n",
    "\n",
    "If you want more technical detail than I present here, you can read my\n",
    "paper on this work at <http://arxiv.org/abs/1203.6890>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Kidney Tumor problem\n",
    "\n",
    "I am a frequent reader and occasional contributor to the online\n",
    "statistics forum at <http://reddit.com/r/statistics>. In November 2011,\n",
    "I read the following message:\n",
    "\n",
    "> “I have Stage IV Kidney Cancer and am trying to determine if the\n",
    "> cancer formed before I retired from the military. ... Given the dates\n",
    "> of retirement and detection is it possible to determine when there was\n",
    "> a 50/50 chance that I developed the disease? Is it possible to\n",
    "> determine the probability on the retirement date? My tumor was 15.5 cm\n",
    "> x 15 cm at detection. Grade II.”\n",
    "\n",
    "I contacted the author of the message and got more information; I\n",
    "learned that veterans get different benefits if it is “more likely than\n",
    "not” that a tumor formed while they were in military service (among\n",
    "other considerations).\n",
    "\n",
    "Because renal tumors grow slowly, and often do not cause symptoms, they\n",
    "are sometimes left untreated. As a result, doctors can observe the rate\n",
    "of growth for untreated tumors by comparing scans from the same patient\n",
    "at different times. Several papers have reported these growth rates.\n",
    "\n",
    "I collected data from a paper by Zhang et al[^3]. I contacted the\n",
    "authors to see if I could get raw data, but they refused on grounds of\n",
    "medical privacy. Nevertheless, I was able to extract the data I needed\n",
    "by printing one of their graphs and measuring it with a ruler.\n",
    "\n",
    "![CDF of RDT in doublings per year.](figs/kidney2.pdf)\n",
    "\n",
    "[fig.kidney2]\n",
    "\n",
    "They report growth rates in reciprocal doubling time (RDT), which is in\n",
    "units of doublings per year. So a tumor with $RDT=1$ doubles in volume\n",
    "each year; with $RDT=2$ it quadruples in the same time, and with\n",
    "$RDT=-1$, it halves. Figure [fig.kidney2] shows the distribution of RDT\n",
    "for 53 patients.\n",
    "\n",
    "The squares are the data points from the paper; the line is a model I\n",
    "fit to the data. The positive tail fits an exponential distribution\n",
    "well, so I used a mixture of two exponentials."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A simple model\n",
    "\n",
    "It is usually a good idea to start with a simple model before trying\n",
    "something more challenging. Sometimes the simple model is sufficient for\n",
    "the problem at hand, and if not, you can use it to validate the more\n",
    "complex model.\n",
    "\n",
    "For my simple model, I assume that tumors grow with a constant doubling\n",
    "time, and that they are three-dimensional in the sense that if the\n",
    "maximum linear measurement doubles, the volume is multiplied by eight.\n",
    "\n",
    "I learned from my correspondent that the time between his discharge from\n",
    "the military and his diagnosis was 3291 days (about 9 years). So my\n",
    "first calculation was, “If this tumor grew at the median rate, how big\n",
    "would it have been at the date of discharge?”\n",
    "\n",
    "The median volume doubling time reported by Zhang et al is 811 days.\n",
    "Assuming 3-dimensional geometry, the doubling time for a linear measure\n",
    "is three times longer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# time between discharge and diagnosis, in days \n",
    "interval = 3291.0\n",
    "\n",
    "# doubling time in linear measure is doubling time in volume * 3\n",
    "dt = 811.0 * 3\n",
    "\n",
    "# number of doublings since discharge\n",
    "doublings = interval / dt\n",
    "\n",
    "# how big was the tumor at time of discharge (diameter in cm)\n",
    "d1 = 15.5\n",
    "d0 = d1 / 2.0 ** doublings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can download the code in this chapter from\n",
    "<http://thinkbayes.com/kidney.py>. For more information see\n",
    "Section [download].\n",
    "\n",
    "The result, `d0`, is about 6 cm. So if this tumor formed\n",
    "after the date of discharge, it must have grown substantially faster\n",
    "than the median rate. Therefore I concluded that it is “more likely than\n",
    "not” that this tumor formed before the date of discharge.\n",
    "\n",
    "In addition, I computed the growth rate that would be implied if this\n",
    "tumor had formed after the date of discharge. If we assume an initial\n",
    "size of 0.1 cm, we can compute the number of doublings to get to a final\n",
    "size of 15.5 cm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "150.76707583570433 2.4209529698496777\n"
     ]
    }
   ],
   "source": [
    "from math import log2\n",
    "\n",
    "# assume an initial linear measure of 0.1 cm\n",
    "d0 = 0.1\n",
    "d1 = 15.5\n",
    "\n",
    "# how many doublings would it take to get from d0 to d1\n",
    "doublings = log2(d1 / d0)\n",
    "\n",
    "# what linear doubling time does that imply?\n",
    "dt = interval / doublings\n",
    "\n",
    "# compute the volumetric doubling time and RDT\n",
    "vdt = dt / 3\n",
    "rdt = 365 / vdt\n",
    "print(vdt, rdt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`dt` is linear doubling time, so `vdt` is\n",
    "volumetric doubling time, and `rdt` is reciprocal doubling\n",
    "time.\n",
    "\n",
    "The number of doublings, in linear measure, is 7.3, which implies an RDT\n",
    "of 2.4. In the data from Zhang et al, only 20% of tumors grew this fast\n",
    "during a period of observation. So again, I concluded that is “more\n",
    "likely than not” that the tumor formed prior to the date of discharge.\n",
    "\n",
    "These calculations are sufficient to answer the question as posed, and\n",
    "on behalf of my correspondent, I wrote a letter explaining my\n",
    "conclusions to the Veterans’ Benefit Administration.\n",
    "\n",
    "Later I told a friend, who is an oncologist, about my results. He was\n",
    "surprised by the growth rates observed by Zhang et al, and by what they\n",
    "imply about the ages of these tumors. He suggested that the results\n",
    "might be interesting to researchers and doctors.\n",
    "\n",
    "But in order to make them useful, I wanted a more general model of the\n",
    "relationship between age and size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A more general model\n",
    "\n",
    "Given the size of a tumor at time of diagnosis, it would be most useful\n",
    "to know the probability that the tumor formed before any given date; in\n",
    "other words, the distribution of ages.\n",
    "\n",
    "To find it, I run simulations of tumor growth to get the distribution of\n",
    "size conditioned on age. Then we can use a Bayesian approach to get the\n",
    "distribution of age conditioned on size.\n",
    "\n",
    "The simulation starts with a small tumor and runs these steps:\n",
    "\n",
    "1.  Choose a growth rate from the distribution of RDT.\n",
    "\n",
    "2.  Compute the size of the tumor at the end of an interval.\n",
    "\n",
    "3.  Record the size of the tumor at each interval.\n",
    "\n",
    "4.  Repeat until the tumor exceeds the maximum relevant size.\n",
    "\n",
    "For the initial size I chose 0.3 cm, because carcinomas smaller than\n",
    "that are less likely to be invasive and less likely to have the blood\n",
    "supply needed for rapid growth (see\n",
    "<http://en.wikipedia.org/wiki/Carcinoma_in_situ>).\n",
    "\n",
    "I chose an interval of 245 days (about 8 months) because that is the\n",
    "median time between measurements in the data source.\n",
    "\n",
    "For the maximum size I chose 20 cm. In the data source, the range of\n",
    "observed sizes is 1.0 to 12.0 cm, so we are extrapolating beyond the\n",
    "observed range at each end, but not by far, and not in a way likely to\n",
    "have a strong effect on the results.\n",
    "\n",
    "![Simulations of tumor growth, size vs. time.](figs/kidney4.pdf)\n",
    "\n",
    "[fig.kidney4]\n",
    "\n",
    "The simulation is based on one big simplification: the growth rate is\n",
    "chosen independently during each interval, so it does not depend on age,\n",
    "size, or growth rate during previous intervals.\n",
    "\n",
    "In Section [serial] I review these assumptions and consider more\n",
    "detailed models. But first let’s look at some examples.\n",
    "\n",
    "Figure [fig.kidney4] shows the size of simulated tumors as a function of\n",
    "age. The dashed line at 10 cm shows the range of ages for tumors at that\n",
    "size: the fastest-growing tumor gets there in 8 years; the slowest takes\n",
    "more than 35.\n",
    "\n",
    "I am presenting results in terms of linear measurements, but the\n",
    "calculations are in terms of volume. To convert from one to the other,\n",
    "again, I use the volume of a sphere with the given diameter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "Here is the kernel of the simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from kidney import Volume\n",
    "\n",
    "def MakeSequence(rdt_seq, v0=0.01, interval=0.67, vmax=Volume(20.0)):\n",
    "    seq = v0,\n",
    "    age = 0\n",
    "\n",
    "    for rdt in rdt_seq:\n",
    "        age += interval\n",
    "        final, seq = ExtendSequence(age, seq, rdt, interval)\n",
    "        if final > vmax:\n",
    "            break\n",
    "\n",
    "    return seq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`rdt_seq` is an iterator that yields random values from the CDF of\n",
    "growth rate. `v0` is the initial volume in mL.\n",
    "`interval` is the time step in years. `vmax` is\n",
    "the final volume corresponding to a linear measurement of 20 cm.\n",
    "\n",
    "`Volume` converts from linear measurement in cm to volume in\n",
    "mL, based on the simplification that the tumor is a sphere:\n",
    "\n",
    "    def Volume(diameter, factor=4*math.pi/3):\n",
    "        return factor * (diameter/2.0)**3\n",
    "\n",
    "`ExtendSequence` computes the volume of the tumor at the end\n",
    "of the interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ExtendSequence(age, seq, rdt, interval):\n",
    "    initial = seq[-1]\n",
    "    doublings = rdt * interval\n",
    "    final = initial * 2**doublings\n",
    "    new_seq = seq + (final,)\n",
    "    cache.Add(age, new_seq, rdt)\n",
    "\n",
    "    return final, new_seq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`age` is the age of the tumor at the end of the interval.\n",
    "`seq` is a tuple that contains the volumes so far.\n",
    "`rdt` is the growth rate during the interval, in doublings\n",
    "per year. `interval` is the size of the time step in years.\n",
    "\n",
    "The return values are `final`, the volume of the tumor at the\n",
    "end of the interval, and `new_seq`, a new tuple containing the volumes\n",
    "in `seq` plus the new volume `final`.\n",
    "\n",
    "`Cache.Add` records the age and size of each tumor at the end\n",
    "of each interval, as explained in the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Caching the joint distribution\n",
    "\n",
    "![Joint distribution of age and tumor size.](figs/kidney8.pdf)\n",
    "\n",
    "[fig.kidney8]\n",
    "\n",
    "Here’s how the cache works.\n",
    "\n",
    "```python\n",
    "class Cache(object):\n",
    "\n",
    "    def __init__(self):\n",
    "        self.joint = thinkbayes.Joint()\n",
    "```\n",
    "\n",
    "`joint` is a joint Pmf that records the frequency of each\n",
    "age-size pair, so it approximates the joint distribution of age and\n",
    "size.\n",
    "\n",
    "At the end of each simulated interval, `ExtendSequence` calls\n",
    "`Add`:\n",
    "\n",
    "```python\n",
    "    # class Cache\n",
    "\n",
    "        def Add(self, age, seq):\n",
    "            final = seq[-1]\n",
    "            cm = Diameter(final)\n",
    "            bucket = round(CmToBucket(cm))\n",
    "            self.joint.Incr((age, bucket))\n",
    "```\n",
    "\n",
    "Again, `age` is the age of the tumor, and `seq` is\n",
    "the sequence of volumes so far.\n",
    "\n",
    "![Distributions of age, conditioned on size.](figs/kidney6.pdf)\n",
    "\n",
    "[fig.kidney6]\n",
    "\n",
    "Before adding the new data to the joint distribution, we use\n",
    "`Diameter` to convert from volume to diameter in centimeters:\n",
    "\n",
    "    def Diameter(volume, factor=3/math.pi/4, exp=1/3.0):\n",
    "        return 2 * (factor * volume) ** exp\n",
    "\n",
    "And `CmToBucket` to convert from centimeters to a discrete\n",
    "bucket number:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "def CmToBucket(x, factor=10):\n",
    "    return factor * math.log(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The buckets are equally spaced on a log scale. Using\n",
    "`factor=10` yields a reasonable number of buckets; for\n",
    "example, 1 cm maps to bucket 0 and 10 cm maps to bucket 23.\n",
    "\n",
    "After running the simulations, we can plot the joint distribution as a\n",
    "pseudocolor plot, where each cell represents the number of tumors\n",
    "observed at a given size-age pair. Figure [fig.kidney8] shows the joint\n",
    "distribution after 1000 simulations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conditional distributions\n",
    "\n",
    "![Percentiles of tumor age as a function of size.](figs/kidney7.pdf)\n",
    "\n",
    "[fig.kidney7]\n",
    "\n",
    "By taking a vertical slice from the joint distribution, we can get the\n",
    "distribution of sizes for any given age. By taking a horizontal slice,\n",
    "we can get the distribution of ages conditioned on size.\n",
    "\n",
    "Here’s the code that reads the joint distribution and builds the\n",
    "conditional distribution for a given size.\n",
    "\n",
    "```python\n",
    "# class Cache\n",
    "\n",
    "    def ConditionalCdf(self, bucket):\n",
    "        pmf = self.joint.Conditional(0, 1, bucket)\n",
    "        cdf = pmf.MakeCdf()\n",
    "        return cdf\n",
    "```\n",
    "\n",
    "`bucket` is the integer bucket number corresponding to tumor size.\n",
    "`Joint.Conditional` computes the PMF of age conditioned on\n",
    "`bucket`. The result is the CDF of age conditioned on\n",
    "`bucket`.\n",
    "\n",
    "Figure [fig.kidney6] shows several of these CDFs, for a range of sizes.\n",
    "To summarize these distributions, we can compute percentiles as a\n",
    "function of size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from kidney import Cache\n",
    "\n",
    "percentiles = [95, 75, 50, 25, 5]\n",
    "cache = Cache()\n",
    "\n",
    "for bucket in cache.GetBuckets():\n",
    "    cdf = ConditionalCdf(bucket)      \n",
    "    ps = [cdf.Percentile(p) for p in percentiles]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Figure [fig.kidney7] shows these percentiles for each size bucket. The\n",
    "data points are computed from the estimated joint distribution. In the\n",
    "model, size and time are discrete, which contributes numerical errors,\n",
    "so I also show a least squares fit for each sequence of percentiles."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Serial Correlation\n",
    "\n",
    "The results so far are based on a number of modeling decisions; let’s\n",
    "review them and consider which ones are the most likely sources of\n",
    "error:\n",
    "\n",
    "-   To convert from linear measure to volume, we assume that tumors are\n",
    "    approximately spherical. This assumption is probably fine for tumors\n",
    "    up to a few centimeters, but not for very large tumors.\n",
    "\n",
    "-   The distribution of growth rates in the simulations are based on a\n",
    "    continuous model we chose to fit the data reported by Zhang et al,\n",
    "    which is based on 53 patients. The fit is only approximate and, more\n",
    "    importantly, a larger sample would yield a different distribution.\n",
    "\n",
    "-   The growth model does not take into account tumor subtype or grade;\n",
    "    this assumption is consistent with the conclusion of Zhang et al:\n",
    "    “Growth rates in renal tumors of different sizes, subtypes and\n",
    "    grades represent a wide range and overlap substantially.” But with a\n",
    "    larger sample, a difference might become apparent.\n",
    "\n",
    "-   The distribution of growth rate does not depend on the size of the\n",
    "    tumor. This assumption would not be realistic for very small and\n",
    "    very large tumors, whose growth is limited by blood supply.\n",
    "\n",
    "    But tumors observed by Zhang et al ranged from 1 to 12 cm, and they\n",
    "    found no statistically significant relationship between size and\n",
    "    growth rate. So if there is a relationship, it is likely to be weak,\n",
    "    at least in this size range.\n",
    "\n",
    "-   In the simulations, growth rate during each interval is independent\n",
    "    of previous growth rates. In reality it is plausible that tumors\n",
    "    that have grown quickly in the past are more likely to grow quickly.\n",
    "    In other words, there is probably a serial correlation in growth\n",
    "    rate.\n",
    "\n",
    "Of these, the first and last seem the most problematic. I’ll investigate\n",
    "serial correlation first, then come back to spherical geometry.\n",
    "\n",
    "To simulate correlated growth, I wrote a generator[^4] that yields a\n",
    "correlated series from a given Cdf. Here’s how the algorithm works:\n",
    "\n",
    "1.  Generate correlated values from a Gaussian distribution. This is\n",
    "    easy to do because we can compute the distribution of the next value\n",
    "    conditioned on the previous value.\n",
    "\n",
    "2.  Transform each value to its cumulative probability using the\n",
    "    Gaussian CDF.\n",
    "\n",
    "3.  Transform each cumulative probability to the corresponding value\n",
    "    using the given Cdf.\n",
    "\n",
    "Here’s what that looks like in code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def CorrelatedGenerator(cdf, rho):\n",
    "    x = random.gauss(0, 1)\n",
    "    yield Transform(x)\n",
    "\n",
    "    sigma = math.sqrt(1 - rho**2);    \n",
    "    while True:\n",
    "        x = random.gauss(x * rho, sigma)\n",
    "        yield Transform(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`cdf` is the desired Cdf; `rho` is the desired\n",
    "correlation. The values of `x` are Gaussian;\n",
    "`Transform` converts them to the desired distribution.\n",
    "\n",
    "The first value of `x` is Gaussian with mean 0 and standard\n",
    "deviation 1. For subsequent values, the mean and standard deviation\n",
    "depend on the previous value. Given the previous `x`, the\n",
    "mean of the next value is `x \\* rho`, and the variance is\n",
    "`1 - rho\\*\\*2`.\n",
    "\n",
    "`Transform` maps from each Gaussian value, `x`, to\n",
    "a value from the given Cdf, `y`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Transform(x):\n",
    "    p = thinkbayes.GaussianCdf(x)\n",
    "    y = cdf.Value(p)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`GaussianCdf` computes the CDF of the standard Gaussian\n",
    "distribution at `x`, returning a cumulative probability.\n",
    "`Cdf.Value` maps from a cumulative probability to the\n",
    "corresponding value in `cdf`.\n",
    "\n",
    "Depending on the shape of `cdf`, information can be lost in\n",
    "transformation, so the actual correlation might be lower than\n",
    "`rho`. For example, when I generate 10000 values from the\n",
    "distribution of growth rates with `rho=0.4`, the actual\n",
    "correlation is 0.37. But since we are guessing at the right correlation\n",
    "anyway, that’s close enough.\n",
    "\n",
    "Remember that `MakeSequence` takes an iterator as an\n",
    "argument. That interface allows it to work with different generators:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'cdf' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-10-3da5dc3eedc2>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;32mfrom\u001b[0m \u001b[0mkidney\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mUncorrelatedGenerator\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mCorrelatedGenerator\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m \u001b[0miterator\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mUncorrelatedGenerator\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcdf\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      4\u001b[0m \u001b[0mseq1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mMakeSequence\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0miterator\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'cdf' is not defined"
     ]
    }
   ],
   "source": [
    "from kidney import UncorrelatedGenerator, CorrelatedGenerator\n",
    "\n",
    "iterator = UncorrelatedGenerator(cdf)\n",
    "seq1 = MakeSequence(iterator)\n",
    "\n",
    "iterator = CorrelatedGenerator(cdf, rho)\n",
    "seq2 = MakeSequence(iterator)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, `seq1` and `seq2` are drawn from\n",
    "the same distribution, but the values in `seq1` are\n",
    "uncorrelated and the values in `seq2` are correlated with a\n",
    "coefficient of approximately `rho`.\n",
    "\n",
    "Now we can see what effect serial correlation has on the results; the\n",
    "following table shows percentiles of age for a 6 cm tumor, using the\n",
    "uncorrelated generator and a correlated generator with target\n",
    "$\\rho = 0.4$.\n",
    "\n",
    "Correlation makes the fastest growing tumors faster and the slowest\n",
    "slower, so the range of ages is wider. The difference is modest for low\n",
    "percentiles, but for the 95th percentile it is more than 6 years. To\n",
    "compute these percentiles precisely, we would need a better estimate of\n",
    "the actual serial correlation.\n",
    "\n",
    "However, this model is sufficient to answer the question we started\n",
    "with: given a tumor with a linear dimension of 15.5 cm, what is the\n",
    "probability that it formed more than 8 years ago?\n",
    "\n",
    "Here’s the code:\n",
    "\n",
    "```python\n",
    "# class Cache\n",
    "\n",
    "    def ProbOlder(self, cm, age):\n",
    "        bucket = CmToBucket(cm)\n",
    "        cdf = self.ConditionalCdf(bucket)\n",
    "        p = cdf.Prob(age)\n",
    "        return 1-p\n",
    "```\n",
    "\n",
    "`cm` is the size of the tumor; `age` is the age\n",
    "threshold in years. `ProbOlder` converts size to a bucket\n",
    "number, gets the Cdf of age conditioned on bucket, and computes the\n",
    "probability that age exceeds the given value.\n",
    "\n",
    "With no serial correlation, the probability that a 15.5 cm tumor is\n",
    "older than 8 years is 0.999, or almost certain. With correlation 0.4,\n",
    "faster-growing tumors are more likely, but the probability is still\n",
    "0.995. Even with correlation 0.8, the probability is 0.978.\n",
    "\n",
    "Another likely source of error is the assumption that tumors are\n",
    "approximately spherical. For a tumor with linear dimensions 15.5 x 15\n",
    "cm, this assumption is probably not valid. If, as seems likely, a tumor\n",
    "this size is relatively flat, it might have the same volume as a 6 cm\n",
    "sphere. With this smaller volume and correlation 0.8, the probability of\n",
    "age greater than 8 is still 95%.\n",
    "\n",
    "So even taking into account modeling errors, it is unlikely that such a\n",
    "large tumor could have formed less than 8 years prior to the date of\n",
    "diagnosis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Discussion\n",
    "\n",
    "Well, we got through a whole chapter without using Bayes’s theorem or\n",
    "the `Suite` class that encapsulates Bayesian updates. What\n",
    "happened?\n",
    "\n",
    "One way to think about Bayes’s theorem is as an algorithm for inverting\n",
    "conditional probabilities. Given $\\mathrm{p}(B|A)$, we can\n",
    "compute $\\mathrm{p}(A|B)$, provided we know\n",
    "$\\mathrm{p}(A)$ and $\\mathrm{p}(B)$. Of course\n",
    "this algorithm is only useful if, for some reason, it is easier to\n",
    "compute $\\mathrm{p}(B|A)$ than\n",
    "$\\mathrm{p}(A|B)$.\n",
    "\n",
    "In this example, it is. By running simulations, we can estimate the\n",
    "distribution of size conditioned on age, or\n",
    "$\\mathrm{p}(size|age)$. But it is harder to get the\n",
    "distribution of age conditioned on size, or\n",
    "$\\mathrm{p}(age|size)$. So this seems like a perfect\n",
    "opportunity to use Bayes’s theorem.\n",
    "\n",
    "The reason I didn’t is computational efficiency. To estimate\n",
    "$\\mathrm{p}(size|age)$ for any given size, you have to run\n",
    "a lot of simulations. Along the way, you end up computing\n",
    "$\\mathrm{p}(size|age)$ for a lot of sizes. In fact, you end\n",
    "up computing the entire joint distribution of size and age,\n",
    "$\\mathrm{p}(size, age)$.\n",
    "\n",
    "And once you have the joint distribution, you don’t really need Bayes’s\n",
    "theorem, you can extract $\\mathrm{p}(age|size)$ by taking\n",
    "slices from the joint distribution, as demonstrated in\n",
    "`ConditionalCdf`.\n",
    "\n",
    "So we side-stepped Bayes, but he was with us in spirit."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
